### This script will extract all static spatial data for all sites for BRT analysis
### Script was appended to add dem, slope, and aspect to 'final' dataset

# Load library SDMTools

library('SDMTools')

# Establish directories

in.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/'
tdir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/AWAPASCII/RawAWAP/'
tmaxdir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/AWAPASCII/RawAWAP/Tmax/'
tmindir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/AWAPASCII/RawAWAP/Tmin/'
rrdir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/Realized_Radiation/'
insoldir = '/home1/99/jc152199/GRASS/ASCIIoutput/'
setwd(in.dir)

# Read in XY data and microclimate data

xy = read.csv('xy.csv',header=T)

# Recreate xy from dataframe final

xy2 = data.frame(unique(cbind(as.character(final$site),as.numeric(final$lat), as.numeric(final$long)))) # Unique xy's to use to bind dem, slope, aspect

names(xy2)[1]='site'
names(xy2)[2]='lat'
names(xy2)[3]='long'

xy2$site = as.character(xy2$site)
xy2$lat = as.numeric(as.character(xy2$lat))
xy2$long = as.numeric(as.character(xy2$long))

#Read in data frame with model data

micro = read.csv('micro_min_max.csv', header=T)

micro=micro[-which(substr(micro$date,5,8)=='0229'),] # Remove records from February 29th

final = read.csv('/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/final.csv') # Dataset that already has d2c, d2s, fpc, realrad, and insol


# Read in static ASCII files (D2C, D2S, FPC)

d2c = read.asc('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/Coastdist/ASCII/coastdist_WTplusbuffer_LatLong_WGS1984_250mres.asc')
d2s = read.asc('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/StreamDist/ASCII/logplusonestreamdist_WTplusbuffer_LatLong_WGS1984_250mres.asc')
fpc = read.asc('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/FPC/ASCII/fpc_WTplusbuffer_LatLong_WGS1984_250mres.asc')
dem = read.asc('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/DEM/ASCII/dem_WTplusbuffer_LatLong_WGS1984_250mres.asc')
slope = read.asc('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/Slope/slope_WTplusbuffer_LatLong_WGS1984_250mres.asc')
aspect = read.asc('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/Aspect/aspect_WTplusbuffer_LatLong_WGS1984_250mres.asc')

# Extract data for static ASCII's

xy$d2c = extract.data(cbind(xy$lat,xy$long),d2c)
xy$d2s = extract.data(cbind(xy$lat,xy$long),d2s)
xy$fpc = extract.data(cbind(xy$lat,xy$long),fpc)

xy2$dem = extract.data(cbind(xy2$lat,xy2$long),dem)
xy2$slope = extract.data(cbind(xy2$lat,xy2$long),slope)
xy2$aspect = extract.data(cbind(xy2$lat,xy2$long),aspect)


# Merge xy with micro

micro_xy = merge(micro,xy,by=c('site'))

# Merge xy2 with final

micro_xy2 = merge(final,xy2,by=c('site'))

# Remove unnecessary columns from micro_xy2, and rename columns appropriately, then write out as a .csv file

final2 = data.frame(micro_xy2[,c(1:13,16:18)])
names(final2)[5]='lat'
names(final2)[6]='long'

write.csv(x=final2, file='/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/final_new.csv', row.names=F)

# Now need a loop to extract realized radiation, insolation, tmin, and tmax values

list.insol = list.files(insoldir, pattern='insol', full.names=T) # List of insolation files
DOY.list = sort(unique(substr(micro$date,5,8))) # List of DOY's format ('mmdd')

micro_xy$insol = NA # Make a blank column for insolation data

for (DOY in DOY.list)

	{
	
	t.pos = which(DOY.list==DOY) # Get the DOY (0-365) from the 'mmdd' DOY 
	
	t.insol = read.asc(paste(list.insol[t.pos],sep='')) # Read in the appropriate insolation ASCII
	
	micro_xy$insol[which(as.character(substr(micro_xy$date,5,8))==DOY)]=extract.data(cbind(micro_xy[which(as.character(substr(micro_xy$date,5,8))==DOY),c(5:6)]),t.insol) # Append insolation data to micro_xy
	
	cat(t.pos,'\n') # Display progress of loop
	
	}
	
# End Loop

write.csv(micro_xy,'micro_min_max_xy_static_insol.csv',row.names=F)












